Simulations of complementarity systems using time-stepping

One of the useful outcomes of the theory of complementarity systems is a new family of methods for numerical simulation of discontinuous systems. Here we will demonstrate the essence by introducing the method of time-stepping. And we do it by means of an example.

Example 1 (Simulation using time-stepping) Consider the following discontinuous dynamical system in \mathbb R^2: \begin{aligned} \dot x_1 &= -\operatorname{sign} x_1 + 2 \operatorname{sign} x_2\\ \dot x_2 &= -2\operatorname{sign} x_1 -\operatorname{sign} x_2. \end{aligned}

The state portrait is in Fig. 1.

f₁(x) = -sign(x[1]) + 2*sign(x[2])
f₂(x) = -2*sign(x[1]) - sign(x[2])
f(x) = [f₁(x), f₂(x)]

using CairoMakie
fig = Figure(size = (800, 800),fontsize=20)
ax = Axis(fig[1, 1], xlabel = "x₁", ylabel = "x₂")
streamplot!(ax,(x₁,x₂)->Point2f(f([x₁,x₂])), -2.0..2.0, -2.0..2.0, colormap = :magma)
fig

Precompiling CairoMakie...
    612.6 ms  ✓ OffsetArrays
    507.1 ms  ✓ OffsetArrays → OffsetArraysAdaptExt
    459.7 ms  ✓ StackViews
    529.4 ms  ✓ PaddedViews
    592.3 ms  ✓ MosaicViews
   1433.7 ms  ✓ Interpolations
    907.9 ms  ✓ Interpolations → InterpolationsUnitfulExt
   1079.0 ms  ✓ KernelDensity
   8627.1 ms  ✓ ImageCore
   1617.2 ms  ✓ ImageBase
   2200.1 ms  ✓ Sixel
   2442.6 ms  ✓ JpegTurbo
   2783.7 ms  ✓ PNGFiles
   1482.5 ms  ✓ WebP
   1199.4 ms  ✓ ImageAxes
    778.6 ms  ✓ ImageMetadata
   1152.1 ms  ✓ Netpbm
  27596.9 ms  ✓ TiffImages
  65816.2 ms  ✓ Makie
  23175.5 ms  ✓ CairoMakie
  20 dependencies successfully precompiled in 118 seconds. 255 already precompiled.

Figure 1: State portrait of the discontinuous system

One particular (vector) state trajectory is in Fig. 2.

using DifferentialEquations

function f!(dx,x,p,t)
    dx[1] = -sign(x[1]) + 2*sign(x[2])
    dx[2] = -2*sign(x[1]) - sign(x[2])
end

x0 = [-1.0, 1.0]
tfinal = 2.0
tspan = (0.0,tfinal)
prob = ODEProblem(f!,x0,tspan)
sol = solve(prob)

using Plots
Plots.plot(sol,xlabel="t",ylabel="x",label=false,lw=3)

Precompiling DifferentialEquations...
    320.8 ms  ✓ StaticArrayInterface → StaticArrayInterfaceOffsetArraysExt
    636.8 ms  ✓ Sparspak
  12044.4 ms  ✓ LoopVectorization
    781.0 ms  ✓ LoopVectorization → SpecialFunctionsExt
    855.3 ms  ✓ LoopVectorization → ForwardDiffExt
    853.0 ms  ✓ TriangularSolve
   4506.7 ms  ✓ RecursiveFactorization
  13509.2 ms  ✓ LinearSolve
   2045.7 ms  ✓ LinearSolve → LinearSolveRecursiveArrayToolsExt
   2075.7 ms  ✓ LinearSolve → LinearSolveBandedMatricesExt
   2089.7 ms  ✓ LinearSolve → LinearSolveEnzymeExt
   2161.6 ms  ✓ NonlinearSolveBase → NonlinearSolveBaseLinearSolveExt
   2169.1 ms  ✓ LinearSolve → LinearSolveFastAlmostBandedMatricesExt
   2523.2 ms  ✓ OrdinaryDiffEqDifferentiation
   3907.6 ms  ✓ OrdinaryDiffEqExtrapolation
  11369.7 ms  ✓ NonlinearSolveFirstOrder
  11997.4 ms  ✓ OrdinaryDiffEqRosenbrock
   3765.9 ms  ✓ BoundaryValueDiffEqCore
  25735.2 ms  ✓ NonlinearSolve
   3877.0 ms  ✓ NonlinearSolve → NonlinearSolveBandedMatricesExt
   3878.6 ms  ✓ NonlinearSolve → NonlinearSolveNLsolveExt
   4730.2 ms  ✓ OrdinaryDiffEqNonlinearSolve
   3859.6 ms  ✓ OrdinaryDiffEqStabilizedIRK
   4120.6 ms  ✓ OrdinaryDiffEqPDIRK
   4835.3 ms  ✓ OrdinaryDiffEqSDIRK
   3541.2 ms  ✓ OrdinaryDiffEqIMEXMultistep
   3961.9 ms  ✓ OrdinaryDiffEqFIRK
   3302.9 ms  ✓ OrdinaryDiffEqExponentialRK
   8459.0 ms  ✓ OrdinaryDiffEqBDF
  45279.8 ms  ✓ BoundaryValueDiffEqMIRK
  16522.8 ms  ✓ OrdinaryDiffEqDefault
   5107.7 ms  ✓ OrdinaryDiffEq
  56868.6 ms  ✓ BoundaryValueDiffEqFIRK
   4447.0 ms  ✓ DelayDiffEq
   5782.5 ms  ✓ BoundaryValueDiffEqShooting
   6908.4 ms  ✓ StochasticDiffEq
   8007.4 ms  ✓ BoundaryValueDiffEq
   6748.5 ms  ✓ DifferentialEquations
  38 dependencies successfully precompiled in 123 seconds. 278 already precompiled.
Precompiling SciMLBaseMakieExt...
   4202.7 ms  ✓ SciMLBase → SciMLBaseMakieExt
  1 dependency successfully precompiled in 5 seconds. 305 already precompiled.

Figure 2: Trajectory of the discontinuous system

We can also plot the trajectory in the state space, as in Fig. 3.

Plots.plot(sol[1,:],sol[2,:],xlabel="x₁",ylabel="x₂",label=false,aspect_ratio=:equal,lw=3,xlims=(-1.2,0.5))

Figure 3: Trajectory of the discontinuous system in the state space

Now, how fast does the solution approach the origin?

Let’s use the 1-norm \|\bm x\|_1 = |x_1| + |x_2| to measure how far the trajectory is from the origin. We then ask: \frac{\mathrm d}{\mathrm dt}\|\bm x\|_1 = ?

We avoid the troubles with nonsmoothness of the absolute value by consider each quadrant separately. Let’s start in the first (upper right) quadrant, that is, x_1>0 and x_2>0, and therefore |x_1| = x_1, \;|x_2| = x_2, and therefore \frac{\mathrm d}{\mathrm dt}\|\bm x\|_1 = \dot x_1 + \dot x_2 = 1 - 3 = -2.

The situation is identical in the other quadrants. And, of course, undefined on the axes.

The conclusion is that the trajectory will hit the origin in finite time: for, say, x_1(0) = 1 and x_2(0) = 1 , the trajectory hits the origin at t=(|x_1(0)|+|x_2(0)|)/2 = 1. But with an infinite number of revolutions around the origin…

How will a standard algoritm for numerical simulation handle this? Let’s have a look at that.

Forward Euler with fixed step size

\begin{aligned} {\color{blue}x_{1,k+1}} &= x_{1,k} + h (-\operatorname{sign} x_{1,k} + 2 \operatorname{sign} x_{2,k})\\ {\color{blue}x_{2,k+1}} &= x_{1,k} + h (-2\operatorname{sign} x_{1,k} - \operatorname{sign} x_{2,k}). \end{aligned}

f(x) = [-sign(x[1]) + 2*sign(x[2]), -2*sign(x[1]) - sign(x[2])]

using LinearAlgebra
N = 1000
x₀ = [-1.0, 1.0]    
x = [x₀]
tfinal = norm(x₀,1)/2
tfinal = 5.0
h = tfinal/N 
t = range(0.0, step=h, stop=tfinal)

for i=1:N
    xnext = x[i] + h*f(x[i]) 
    push!(x,xnext)
end

X = [getindex.(x, i) for i in 1:length(x[1])]

Plots.plot(t,X,lw=3,label=false,xlabel="t",ylabel="x")

Backward Euler

\begin{aligned} {\color{blue} x_{1,k+1}} &= x_{1,k} + h (-\operatorname{sign} {\color{blue}x_{1,k+1}} + 2 \operatorname{sign} {\color{blue}x_{2,k+1}})\\ {\color{blue} x_{2,k+1}} &= x_{1,k} + h (-2\operatorname{sign} {\color{blue}x_{1,k+1}} - \operatorname{sign} {\color{blue}x_{2,k+1}}). \end{aligned}

Formulation using LCP

Instead solving the above nonlinear equations with discontinuities, we introduce new variables u_1 and u_2 as the outputs of the \operatorname{sign} functions: \begin{aligned} {\color{blue} x_{1,k+1}} &= x_{1,k} + h (-{\color{blue}u_{1}} + 2 {\color{blue}u_{2}})\\ {\color{blue} x_{2,k+1}} &= x_{1,k} + h (-2{\color{blue}u_{1}} - {\color{blue}u_{2}}). \end{aligned}

But now we have to enforce the relation between \bm u and \bm x_{k+1}. Recall the standard definition of the \operatorname{sign} function is \operatorname{sign}(x) = \begin{cases} 1 & x>0\\ 0 & x=0\\ -1 & x<0, \end{cases} but we change the definition to a set-valued function \begin{cases} \operatorname{sign}(x) = 1 & x>0\\ \operatorname{sign}(x) \in [-1,1] & x=0\\ \operatorname{sign}(x) = -1 & x<0. \end{cases}

Accordingly, we set the relationship between \bm u and \bm x to \begin{cases} u_1 = 1 & x_1>0\\ u_1 \in [-1,1] & x_1=0\\ u_1 = -1 & x_1<0, \end{cases} and \begin{cases} u_2 = 1 & x_2>0\\ u_2 \in [-1,1] & x_2=0\\ u_2 = -1 & x_2<0. \end{cases}

But these are mixed complementarity contraints we have defined previously! \boxed{ \begin{aligned} \begin{bmatrix} {\color{blue} x_{1,k+1}}\\ {\color{blue} x_{1,k+1}} \end{bmatrix} &= \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} + h \begin{bmatrix} -1 & 2 \\ -2 & -1 \end{bmatrix} \begin{bmatrix} {\color{blue}u_{1}}\\ {\color{blue}u_{2}} \end{bmatrix}\\ -1 \leq {\color{blue} u_1} \leq 1 \quad &\bot \quad -{\color{blue}x_{1,k+1}}\\ -1 \leq {\color{blue} u_2} \leq 1 \quad &\bot \quad -{\color{blue}x_{2,k+1}}. \end{aligned} }

9 possible combinations

There are now 9 possible combinations of the values of u_1 and u_2. Let’s explore some: x_{1,k+1} = x_{2,k+1} = 0, while u_1 \in [-1,1] and u_2 \in [-1,1]:

\begin{aligned} \begin{bmatrix} 0\\ 0 \end{bmatrix} &= \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} + h \begin{bmatrix} -1 & 2 \\ -2 & -1 \end{bmatrix} \begin{bmatrix} {\color{blue}u_{1}}\\ {\color{blue}u_{2}} \end{bmatrix}\\ & -1 \leq {\color{blue} u_1} \leq 1, \quad -1 \leq {\color{blue} u_2} \leq 1 \end{aligned}

How does the set of states from which the next state is zero look like? \begin{aligned} -\begin{bmatrix} -1 & 2 \\ -2 & -1 \end{bmatrix}^{-1} \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} &= h \begin{bmatrix} {\color{blue}u_{1}}\\ {\color{blue}u_{2}} \end{bmatrix}\\ -1 \leq {\color{blue} u_1} \leq 1, \quad -1 &\leq {\color{blue} u_2} \leq 1 \end{aligned}

\begin{bmatrix} -h\\-h \end{bmatrix} \leq \begin{bmatrix} 0.2 & 0.4 \\ -0.4 & 0.2 \end{bmatrix} \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} \leq \begin{bmatrix} h\\ h \end{bmatrix}

For h=0.2

using LazySets
h = 0.2
H1u = HalfSpace([0.2, 0.4], h)
H2u = HalfSpace([-0.4, 0.2], h)
H1l = HalfSpace(-[0.2, 0.4], h)
H2l = HalfSpace(-[-0.4, 0.2], h)

Ha = H1u ∩ H2u ∩ H1l ∩ H2l

using Plots
Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))

Indeed, if the current state is in this rotated square, then the next state will be zero.

Another

u_1 = 1, u_2 = 1:

\begin{aligned} \begin{bmatrix} {\color{blue} x_{1,k+1}}\\ {\color{blue} x_{1,k+1}} \end{bmatrix} &= \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} + h \begin{bmatrix} -1 & 2 \\ -2 & -1 \end{bmatrix} \begin{bmatrix} {1}\\ {1} \end{bmatrix}\\ \color{blue}x_{1,k+1} &\geq 0\\ \color{blue}x_{2,k+1} &\geq 0 \end{aligned} which can be reformatted to \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} + h \begin{bmatrix} -1 & 2 \\ -2 & -1 \end{bmatrix} \begin{bmatrix} 1\\ 1 \end{bmatrix}\geq \bm 0 and further to \begin{bmatrix} x_{1,k}\\ x_{2,k} \end{bmatrix} \geq h \begin{bmatrix} -1\\ 3 \end{bmatrix}

using LazySets
h = 0.2
A = [-1.0 2.0; -2.0 -1.0]
u = [1.0, 1.0]
b = h*A*u

H1 = HalfSpace([-1.0, 0.0], b[1])
H2 = HalfSpace([0.0, -1.0], b[2])
Hb = H1 ∩ H2

using Plots
Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
Plots.plot!(Hb)

All nine regions

using LazySets
h = 0.2
A = [-1.0 2.0; -2.0 -1.0]

u = [1, -1]
b = h*A*u

H1 = HalfSpace(-[1.0, 0.0], b[1])
H2 = HalfSpace([0.0, 1.0], -b[2])
Hc = H1 ∩ H2

u = [-1, 1]
b = h*A*u

H1 = HalfSpace([1.0, 0.0], -b[1])
H2 = HalfSpace(-[0.0, 1.0], b[2])
Hd = H1 ∩ H2

u = [-1, -1]
b = h*A*u

H1 = HalfSpace([1.0, 0.0], -b[1])
H2 = HalfSpace([0.0, 1.0], -b[2])
He = H1 ∩ H2

using Plots
Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
Plots.plot!(Hb)
Plots.plot!(Hc)
Plots.plot!(Hd)
Plots.plot!(He)

Solutions using a MCP solver

M = [-1 2; -2 -1]
h = 2e-1
tfinal = 2.0
N = tfinal/h

x0 = [-1.0, 1.0]
x = [x0]

using JuMP
using PATHSolver

for i = 1:N
    model = Model(PATHSolver.Optimizer)
    set_optimizer_attribute(model, "output", "no")
    set_silent(model)
    @variable(model, -1 <= u[1:2] <= 1)
    @constraint(model, -h*M * u - x[end] ⟂ u)
    optimize!(model)
    push!(x, x[end]+h*M*value.(u))
end

t = range(0.0, step=h, stop=tfinal)
X = [getindex.(x, i) for i in 1:length(x[1])]

using Plots
Plots.plot(Ha, aspect_ratio=:equal,xlabel="x₁",ylabel="x₂",label=false,xlims=(-1.5,1.5),ylims=(-1.5,1.5))
Plots.plot!(Hb)
Plots.plot!(Hc)
Plots.plot!(Hd)
Plots.plot!(He)
Plots.plot!(X[1],X[2],xlabel="x₁",ylabel="x₂",label="Time-stepping",aspect_ratio=:equal,lw=3,markershape=:circle)
Plots.plot!(sol[1,:],sol[2,:],label=false,lw=3)